“Calhoun 


Institutional Archive of the Naval Postgraduate School 





Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations l. Thesis and Dissertation Collection, all items 


1966 


A computer program for evaluating lattice sums. 


Dickmann, David B. 


http://hdl.handle.net/10945/9468 


This publication is a work of the U.S. Government as defined in Title 17, United 
States Code, Section 101. Copyright protection is not available for this work in the 
United States. 


Downloaded from NPS Archive: Calhoun 


| Calhoun is the Naval Postgraduate School's public access digital repository for 
| г D U DLE Y research materials and institutional publications created by the NPS community. 
FW : Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 


"WM KNOX appointed — and published — scholarly author. 


LIBRARY Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 


NPS ARCHIVE 


1966 
DICKMANN, D. 





A COMPUTER PROGRAM FOR 
EVALUATING LATTICE SUMS 


DAVID B. DICKMANN 


Librery 
U. S. Naval Postgraduate Sch. 
Monterey, Calitornia 


- o" EA LIBRARY 
NAVAL POSTGRADUATE SCROOQ! 
MONTEREY CA 93543-5101 E 





This document has been approved for public 
release and:ale: its distribution is uniimited, 








De -—- эь» 








А СО? лг BROGRAMIEBEOR 
EVALUATING LATTICE SUMS 
by 
David B. Dickmann 


Lieutenant, Unitéd States Navy 
B.S., Oregon State University, 1958 


Submitted in partial fulfillment 
for the degree of 
MASTER OF SCIENCE IN PHYSICS 


from the 


UNITED- STATBS NAVAL POSTGRADUATE SCHOOL 


May 1966 


\ 
\ 
А 


ivan, d, ABSTRACT 


A recurring problem in solid state physics is that of 
determining the electrostatic field or potential at a point 
within the crystal lattice arising from the molecules of the 
lattice themselves. A typical problem is that of finding 
the electrostatic binding potential of an ionic lattice, 
assuming point changes. 

Nijboer, de Wette and others have worked toward 
perfecting a method for calculating the electric field 
arising from a lattice of dipoles. This work culminated 
in a paper by de Wette and Schacher detailing a procedure 
which ensures a rapidly converging Summation of the three 
dimensional lattice sums involved. 

The present paper is an extension of the latter, 
presenting a computer program for evaluating dipole fields 
in a crystal lattice using the equations obtained by de Wette 


and Schacher. 
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BACKGROUND 
A recurring problem in solid state physics is that of 
determining the electrostatic field or potential at a point 
within the crystal lattice arising from the molecules of the 
lattice themselves. А typical problem is that of finding 
the electrostatic binding potential of an ionic lattice, 
assuming point changes. 


1,2 and others have worked toward 


Nijboer, de Wette 
perfecting a method for calculating the electric field 
arising from a lattice of dipoles. This work culminated 
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in a paper by de Wette and Schacher- detailing a procedure 
which ensures a rapidly converging summation of the three 
dimensional lattice sums involved. 

The present paper is an extension of the latter, 
presenting a computer program for evaluating dipole fields 


in a crystal lattice using the equations obtained by de Wette 


and Schacher. 


1 
B.R.A. Nijboer and F'.W. de Wette, Physica 1A, 422 


(19583 


NA de Wette, Phys. Rev. 123, 103 (1961) 


I ow. de Wette and G.E. Schacher, Phys. Rev. 137, 
A78 (1965) 


THEORY 
Virtually all problems concerning crystal self fields 
reduce ultimately to evaluating a lattice sum of the 


general form: 


у) = df GCc- ag, V. - Das di CAs) 


where x, y, and z are coordinates of the point at which 
the sum is to be evaluated; f ( ) represents the contri- 
bution to the sum from a particular lattice point; and 
аа), bao, саз аге the coordinates of the lattice points 
whose contribution is being considered (a4, ас, апа аз 
being lattice basis vectors). 

Finding the electric field due to a lattice of dipoles 
reduces to evaluating such a sum, where the summation 
is, mathematically, conditionally convergent. ‘That is, the 
value of the sum depends on the order in which the sum- 
mation is carried out (analagous to the fact that the 
internal field of a crystal depends on the shape of the 
crystal). Thus, any procedure for evaluating dipole sums 
must, at some point, make an assumption of a particular 
crystal shape. 

ST a complete description see (2) 
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In the procedure used here these conditionally 
convergent triple sums over the real lattice are converted 
into absolutely (rapidly) convergent double sums over the 
recipricol lattice. Converting from a triple to a double 
sum forces an assumption of a slab-shaped crystal with 
the short dimension in the z-direction (perpendicular to 
the plane described'by the lattice vectors а, апа а. 

See figure 1 for relations between x, y, z coordinates 
and the lattice basis vectors). 

To achieve the conversion to an absolutely convergent 
double sum, the original sum is first Fourier transformed, 
making use of the property of lattice sums that the sum- 
mation of a function over a lattice in real space is equal 
to the sum of its F'ourier transform over the recipricol 
lattice (in Fourier space). In Fourier space it is then 
possible to carry out the z-summation analytically, leaving 
an absolutely convergent double sum in recipricol space 
(with the assumption of a slab-shaped form). 

A distinction is made between the cases where the 
evaluation point does or does not lie in an x-y plane of 
dipoles (z= сәз Тор 272 саз, с ап integer). In the first 
case the function to be summed has a simple pole (with 
respect to the z-summation) at z-0. In this event an 
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auxiliary convergence function is used and the sum is 
broken up into two sums, one in real space over the x-y 
plane containing the evaluation point and one in recipricol 
Space over the rest of the lattice. Іп the second case 
the function is everywhere regular and is F'ourier trans- 
formed as it stands. 

In either event the method is directly applicable only 
to lattices in which all dipoles are of the same magnitude 
and orientation. ‘To evaluate the field of a lattice contain- 
ing more than one orientation and/or magnitude of dipoles 
it is necessary to divide the lattice into sublattices, each 
containing only one kind of dipole, evaluate the field for 
each sublattice, and combine the fields by Superposition. 

The summation procedure, applied to a sublattice, 
produces nine numbers Su > l,j - х,у,2 $ where Sis is the 
i-component of the electric field (in Coulombs) due to a 
lattice of unit dipoles (one electronic charge and length 
1A) oriented in the j-direction. 

While program inputs are in terms of the normal 
crystallographic quantities (а), а2, аз, ©, 3, 2,58 ) the 
calculations and outputs are in terms of a cartesian 
reference frame related to the crystal basis vectors as 


detailed in figures 1-3. 





Figure 1: Unit Cell: 1 dipole per unit cell volume, 
Basis vectors a4, ар, аҙ. Crystal angles: a< d.,d, ; 
0*4, d, ; "q.d. 5 542,9, = 


x axis parallel to aq, y axis in a] =- ap plane perpendicular 


to x axis, z axis perpendicular to xey plane. 





Figure 2: Point at which field is to be evaluated: 


located with respect to lattice point by fractions of 


basis vectors Q; 7j, d, *j«dz *j,d, ; 0 Sj; */ 


— ————Á— шеф сш EE ње 





Figure 3: 


Decomposition of arbitrary dipole: 


р =р»Т =р, ? +Р К 


— e — A MÀ —À HMM ee cc C —  — че 


Then given a lattice of arbitrary dipoles p, as in figure 


3, we have: 


р, = реше сов @ 
ру = р sin © sin Ø 
Pam ae cos 9 


so that the field at some point (14, ір, j3) is given by: 
Eal) = p,S,, (i) +pySyy (i) tpzS x26) 
By(j) = Am 


The nine Su are not all independent of course. They 


satisfy 
Sij = Sji 
S xx E Syy Е 522 = 
for all lattices. Other special relations (such as S., = 0, 


1J 


i # j, for cubic lattices) are also obtained for lattices 
with a high order of symmetry. 

In a lattice containing more than one kind of dipole 
two special considerations arise in decomposing the lattice 
into sublattices. First and most obviously, one must 
choose the j's in each sublattice, determining the point at 
which the fields are to be evaluated, such that they refer 
to a common point in the composite lattice. The second 
consideration refers back to the shape dependence of the 
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sum. In computing the sums over each sublattice, the 
shape dependence has introduced a depolarization factor 
L = hir/(vol of unit cell in АЗ) into the Sz sum. Then 
to superpose the fields from several sublattices with 
various Етте. e (referred to the composite lattice) 
it is necessary to modify all the Sij to bring all sublattices 
to a common x-y-z orientation by manipulating the 
depolarization factors for each sublattice. A similar 
consideration arises in comparing fields obtained by this 
method and those obtained by some other method which 
assumes a different crystal shape or with experimental 
results on a specific crystal of known shape. Kittel? 
considers the general problem of depolarization factors 
in some detail and tabulates the depolarization factors 
for most commonly used shapes (spherical, cubic, slab, 


etc.) An example of such considerations for a benzene 


model is contained in Appendix D. 


Эс. Kittel, Introduction to Solid State Physics, 
end edition, 158, (John Wiley and Sons, Inc., New York, 


1956) 
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PROGRAM COMPOSITIONÎ 


The actual translation of the method of Schacher 
and de Wétte is relatively straightforward although quite 
lengthy. 

The program is written in the FORTRAN 60 language 
for the CDC 1604 computer and is compatable with most 
FORTRAN monitors. The program is divided into a main 
program, containing all input-output routines and a segment 
to choose a subroutine, and eight subroutines to do the 
actual summation. 

Simplified flow charts of the entire program and 


various segments are shown in figures h-7. 


Appendix A contains the program listing and details 
of timing, storage, etc. 
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Figure 4: Flow Main Program - LATSUM 





nput lattice parameters and convergence 
limits 


5 Па 
Олык aCeee Parameters) 
16 -еггог__ Test lattice parameters for correctness 


нь aaah) 
pt [Compute basic constants 
"Же t 
Examine lattice parameters to determine 
which subroutine to use 


(Subroutine SUM1 
(-—[Subroutine SUMe| 
Do actual | 
summation Subroutine SUM3 


Subroutines SUMĻ4, SUMO, 
GAM1,GAM2 











sSubnoutines| SUM5,SUM6, 
„САМ САМЕ 


„=== +2 ы чл ЧЕ 
V 


L 


Figure 5: Common Flow-Subroutine SUM1, SUM2, SUM3 


Main program uses SUM1 if XJ340, XDEL#0 


SUM2 XJIJI340, XDEL=0,XJ1+XJI2A0 
SUM3 XJ3#0, XDEL=0,XJ1+XJ2=0 


ENTER FRƏM 
LATSUM 


Compute basic constants 












Compute contribution to Sij 
for lattice points on 
А1 апа A2 axis 


Convergence 
tests met? 


Yes 


Compute contribution to 5., 


for lattice points in 
A1-A2 plane 





Convergence 
tests met? 


Yes 


RETURN 
TO LATSUM 


Figure 6: Common Flow - Subroutine SUM4, SUM5, GAMi 


Main program uses SUM) & SUM6 if XJ3=0, XDEL#0 
SUM5 & SUM6 if XJ3=0, XDEL=0 


SUMA or SUMS5 
ENTER BROM 
LATSUM 

» constants 


Compute recip. lattice 
contribution to Si 
for lattice points on 


Ai & A2 axis Subroutine 
GAM1 


Compute auxilliary 
convergence functions 








Convergence 
tests met? 


ГИА Tt x^)? and 
ЈУ СИТА Ја 






Yes 
Compute recip. lattice | 
contribution to Si p uw 
for lattice points in 
A1-A2 plane 





No Convergence 
tests met7 







ГО ЕЛ ЕЛЕГО ТАЕ ТМ 
AND THEN ENTER 


SUBROUTINE SUMO 
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Figure 7: Flow Subroutine SUM6, GAM2 


Used whenever XJ3=0 


SUM6 
PENI Pee evi 
LATSUM 


Compute basic constants 


No "XJ14XJ2-20 
? 


Yes 


Add constant Секта ==: 


Compute real lattice contribution 
to Su for lattice points on A1 Subroutine GAM2 
| ^ | TCompute auxilliary 


convergence functions 
MY2 1X *) апа 
CS mima) 


P 
=> M, 



















Convergence 
Tests met? 






Compute real lattice contribution 
to Su for lattice points in A1—-A2 
| olane 


/ Convergence* 
tests met? 


| Yes 
RETURN ТО 
LATSUM l 
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As can be seen from the flow charts, a deck can be 
made up for a particular lattice sum by appending to the 
main program either subroutines SUM1, SUMe, and SUM3 
or subroutines SUM, SUM5, SUM6, GAMi1 апа САМ? 
depending on whether Ја is not zero ог zero respectively. 
Further, if sums are to be computed for several sets of 
data the main program deck can be modified by looping 
from the input segment to the END card. 

Inputs to the program are given in Table 1. 

Output is scattered throughout the main program. 

A check output of the input parameters always appears. 
If there is an error in the input parameters the following 
is printed and control returned to the monitor: 

FROGRAMNMEBSTOPERPED BY INPUT ÆRROR 

DESIGNATED BY K= 


The individual error codes are: 


K ERROR 

1 А1 0. 

2 А2%0. 

h АЗ <0. 

10 Хот S 1-І 

20 X J2:« 0. о ш 21. 

40 XJ3«0. or XJ32» 1. 

100 XALPH€0. or XALPH 290. 
200 08 Е = 0) (Әс XBET 0. 
4.00 XGAM <0. or XGAM=180. 
1000 XDEL<0. or XDEL> 90. 
2000 N= 0 
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Table 1: 


BABEL 
A1 
AR 


ace 
XJ2 
хХ.23 
XALPH 
хып 
САМ 


X DEL 


AEX 
REX 
А ЕҮ 
REM 
A EZ 
REZ 









Program Inputs 


ЮГ = DESCRIPTION 


Length of basis vector a; 


Firror limits for convergence 
ests. Summation will contin- 
e until the next term to be 
summed is less than 
(x) (x) 
AE(y)+RE(y)times absolute 
(2) 2) 
value of sum at that time 
РОН Ie EI CODntrol e» M ЕЕ 
ABYGREY control 5; уз зу 
AEZEREZ control 5 265, 
fuller discussion is "o Ks 
opendix B 
Summation Cutoff. 





erms in x-y plang if converg- 
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Stops sum 
ation after covering 2N by 2N 






ALLOWED 
EIMITS ; 


ai^ 0 
g £ ji е1 
0 < ХАГ,РН = 90 


2<ХБЕТ900 
)« XGAM « 180 


De XDEL 90 


NONE 





N> 0 





Multiple errors result in a K value which is the sum of the 
individual K-values involved. 

A second error output (and transfer back of control) 
occurs if N is so small that it is impossible to check the 
sums for convergence. This subject is covered more 
completely in appendix B. ‘The printout is self-explanatory: 

PROGRAM STOPPED BECAUSE М TOO SMALL. 

N MUST BE AT LEAST -- + 1 TO ENSURE A 

GOOD TEST FOR CONV ERGENGS, 

These are the only error outputs contained in the 
program other than those that might arise from library 
routines. There are no error stops, pauses, dumps, 
sense switches or sense lights used. 

Normal output (presuming no input errors) includes 
the nine sums Si displayed as a labeled array. Following 
this is one of the statements: 

(SUM1) 

SUMMING WAS DONE BY SUBROUTINE (SUMe2) 

(SUM3) 
(DELTA AND J3 BOTH DIFFERENT 
WHICH ASSUMES FROM ZERO) 
(DBETA=6 BUT J1+J2 AND J3 NOT 
EQUAL TO ZERO) 
(J1, 32, AND DELTA ALL ZERO BUT 
J3 NOT ZERO) 


or 
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(SUM, & SUM6) 
SUMMING WAS DONE BY SUBROUTINES (SUM) & SUM6) 
(SUM5 & SUMÓ) 
(SUMS5 © SUM) 


(23-0 ВОТ DELTA AND J14342 NOT ZERO) 
ASSUMING (431,432, & J3 ALL ZERO BUT DELTA NOT 
ZERO) 
(I3 AND DELTA ZERO BUT J1+J2 NOT 
ZERO) 
(J1; J2, J3 AND DELTA ALL ZERO) 
Following this is a statement of whether or not 


covergence test were met, and if not, a measure of the 


error introduced. 
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E. 


B. 


ER 


APPENDIX A 


PROGRAM LISTINE 


IDENTIFICAT FON 

Title Electric Field Components by de Wette 
- Schacher Equations 

Co-op Id. Z1-NPGS-LATSUM (F-60,63) 

Category All Others 

Programmer David Dickmann 

Date February 1966 

PURPOSE 


This Fortram program computes electric fields per 


unit-dipoles for a lattice of similar dipoles. ‘The 
program accepts lattices of any symmetry including 
triclinic and all higher symmetries, and computes field 
components in a cartesian reference frame, using the 
de Wette-Schacher equations. 


USAGE 


dl 


Calling sequence: ‘This is a main program. See 
input formats. 


Input arguments: 


A1, A2, A3 are the lengths, in Angstroms, of 
the lattice unit cell basis vectors. 


X4J1, XJ2, XJ3 are fractions of the basis 
vectors A1, A2, A3 designating the location of 
the point at which the fields are to be evaluated. 


XALPH is the crystal angle between the A2 and 
A3 axes, in degrees and decimals. 


XBET is the crystal angle between the A1 and 
А.З axes, in degrees and decimals. 
AX QGAM is the crystal angle between the A1 and 
A2 axes, in degrees and decimals. 
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XDEIL is the inclination-of the A3 axis from the 
perpendicular to the A1-A2 plane, in degrees and 
decimals. If unknown, enter as XDEL=90 and 
correct value is generated by program. 


AEX, AEY, AE@ are the absolute error limits 
for the values of the x, y, and z components of 


the electric fields. 


REX, REY, REQ are the relative error limits 
for the corresponding electric field components. 


N is a summation cutoff value. If error limits 

are not met after “~ LN“ terms are summed, the 
program returns the field values at that time 

and a measure of the maximum error introduced 

Бу comEOff. 

Space required: 6446 cells (F-60), 4209 cells (F-63) 
lemporary storage: None 

Alarms: None 


Error returns: 


a. Error in input parameters--Control returned 
to monitor with printetr: 


„АРСИ AM SHOPPED BY INPUT 
ЕС РЕ ЕТ В ВАЙ = 


Individual error codes аге: 


E Error 

1 A120. 

2 А2 2 0. 

д. АЗ<0. 

10 ЕЈ < OL Oron Tl >i, 

20 Ра Је сишог Х<Ја> 1, 

40 ЖЗ НОО ХО > +, 

100 XALPH#0. or XALPH >90. 
200 EE qo Torn XBET >90. 
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4.00 XGAM#0. or XGAM2180. 
1000 XDEL <0. or XDEL>90. 
2000 NO 


Multiple errors result in à K value which is 
the sum of the individual values above. 


b. WN value too small to ensure good test for 
convergence of field sums. Control returned 
to monitor with self-explanatory printout. 


7, Error Stops: None except those associated with 
machine arithmetic errors. 


8. Input and output tape mountings: Normal 
9.. Inpwt and output formats 


a. Input is contained on four cards under the 
following format control: 


READ 1,A1,A2,A3,XJ1,XJ2,XJ3,XALPH, 
XBET,XGAM,XDEL,N,AEX,REX, 
AEY,REY,AEZ,REZ 

1 FORMAT (3F10.0/3F10.0/LE'10.0/I4,,6E10.0) 


b. Output: Various 
(1) Error outputs as per 6 
(2) Su - - the i-th component of the electric 
field due to a unit dipole oriented in 
the j-direction. Printed as an array 


under 1PE19.10 with labels SUMEXX, etc. 


(3) Hollerith prints indicating subroutines 
used and results of convergence tests. 


А у ои 


4" 35 к 
2 , 
? pu у ме м. 


10. Selective stop and jump settings: Not applicable 
1%. Шиле: 


a. Compile time: Дд тіп 20 sec 


p 


b. Run time: Average 1.4 sec per lattice, 
averaged over 336 combinations of lattices, 
error limits, and evaluation points. 


12. Accuracy: Relative error of S; + min; j(REX,REY, 
наты 
Absolute errér of S 

AEZ) 


етіп, (АЕХ,АБҮ, 


1) 
13. Equipment configuration: CDC 1604 with 
FORTRAN compiler or CO-OFP»monitor system 
including F-63 
14,. Library routines used: 
ES CO ISI" 
SINE 
TOSE 
АТАМЕ 
ASINE 
ACOSFE 
EXPE 
ABSF 
FIXFLOAT 
FLOAT FIX 


15. References: 


а. Ш. EV DICEN -—-wmputer Program 


for Evaluating Lattice Sums, Thesis, U.S. 
Naval Postgraduate School, Monterey, 1966 


b. F. W. de Wette and G. E. Schacher, Phys. 
Rev., 137, A78 (1965) 


D. METHOD OF ALGORITHM See references (a) and 
(b). 


E. FLOW CHART See reference (a). 
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АЗЕ БӘМ B 


БОММАТттТОМ СЕБЕБ ЕЈ\ЧС Е ЖАМЫШ ЕЕЕ ШЕ 


The general lattice sums to be evaluated in the 
program are of the form 


oo 2 / -с20/ 
ED "eb h; Cu, ‚ “4; Ø, CU, , s ) $ /-Ф, CM, , Me) ea Uu 


—-00 


я егей; 3 
И, ЛЕСІП” 
/7 6; (U, МА Је ә” 


where C. is a constant, hj (44), ) іѕ a function of £4, 
and A, which generally increases with increasing 4, and 42 
and (f; is a phase factor of the form expt-2fri(c;u, "Ал 
Because of the phase factors it is possible to have a 
distribution of significant terms over the 41-4 plane such 
as shown in figure 8. 

To ensure inclusion of significant terms in such a 
distribution would require either very complex convergence 
tests or inclusion of many more terms then required. 

To avoid this, the sums were converted to semi- 
infinite sums over the first quadrant of the 41-42 plane. 
This led to the possible distributions shown in figure 9, 
where the degree of asymmetry about the 41-42 axis 
depended on the ratio of the lengths of the a, and ас 


basis vectors апа the angle between them, and the 
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A a 
————————  À 2 А”, 
A ^62 


222,222 Hi 


Figure 8: Possible distributions of significant terms of 
sum over,4.-45 plane. Shaded area contains 
all terms of absolute value larger than a given 


constant. 


<< 


"E 


А; 


>A 


M, 


Figure 9: Possible distributions of significant terms of 


converted sums 
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particular sum being evaluated, S; While these distributions 


Г. 
still leave problems .in developing convergence tests, the 
problems are not as formidable as for the original equations. 

The solution used in this program proceeded in 
several steps. First the axial terms, Af,=0, Үр-0, 27 
ee. апа,/--0, 4:1, 2,3 . . . were summed together 
апа the convergence test applied to the sum of each pair, 
2450, Мәті; Ж-і, %--0, until the tests were met for each 
sum Si: or N pairs were summed. 

The non-axial terms were then treated separately. 
First a skew factor, JIG = 2 ж МАХ 22; AI) cos 8 +1, 
was generated. This number is such that the largest term 
of the sum (in absolute value) must occur at some point 
within the box 0#4,4Z242J/IG . This was tested against N 
to ensure N > JIG so that the sums would include the 
largest term. 

Summing was done by summing all terms > = 1, 2, 
ју « « for =l and then тоге == etc. Phe *convergeuce 
test were threefold. First to ensure for given 44 that 
Ada assumed all values up to JIG and that this was done 
for all values of 4, up to JIG. This ensured picking up 
the most significant term of this sum regardless of its 


magnitude. Second to ensure that for given g4,-a if some 
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term, ><а, did not meet the convergence tests then 
summing would be done on the pencil y4,=atl, Mo=1, 2, . ., 
atl. Last, that summing would not exceed the terms 
447N;, 427i; 447); 4427 N, regardless of other tests. 

The flow chart for this testing is shown in figure 10. 
Actual convergence zones for a cubic lattice and a triclinic 
lattice are shown in figures 11 and 12. 

After completing the program it was checked by 
running sums on 112 lattices (distortions of the basic cubic 
lattice) for which sums are known from geometric 
considerations. The error analysis for a relative error 


8 and an absolute 


limit РУХ étc)wof 1074, 107°yeana 107 
error limit (AEX, etc) equal to the relative limit is given 


in Table 2. 
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Figure 10a Flow chart of convergence tests common to 
SUMi-SUM6. Integer variable I refers to 
values of summing index ul, J to 44). 
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Figure 10b Continuation of flow chart 
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"Figure 11: Actual convergence zones plotted for unit 
cubic lattice using library subroutine 
Ө Маг. 
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M, 


d 


Figure 12: Actual convergence zones plotted for a 
triclinic lattice using library subroutine 


CONTOUR; 
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Table 2: 


2а: 


Lattice 
Type 


Triclinic 


Tetragonal 


Trigonal 


Coubic 


Error analysis of known sums S. Lo+-eackrbex the 
left-hand numbers are absolute errors, the right 
hand numbers are relative errors, the top numbers 
are maximum values, the bottom are average values, 
and the center integer is the number of sums 
evaluated. 


Input relative error limit = 1074 
Input absolute error limit = ow 
All errors X 107 


Subroutine Used 
TI SUM2 SUM3 SUML56 SUM56&6 





2b: Input relative error limit = 1076 
Input absolute error limit = 10-0 


All errors Х 1078 





Lattice Subroutine Used 
‘Type SUM1 SUM2 SUM3 сааи SUM566 
Monoclinic 

ee oe — 
Tetragonal Par 





-8 


ес: Input relative error limit = 10 
Input absolute error limit = 1079 
Al errors X 10710 
Lattice Subroutine Used 
‘Type SUM1 SUM2 SUM3 Од SUM556 










Triclinic 


Monoclinic 


Cubic 
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КОБ ГООО (© 


INCOMPLETE GAMMA FUNCTION 


In evaluating the lattice sums when the evaluation 
point lines in the x-y plane an auxilliary convergence 
function is used to smooth and remove the simple pole of 
the function being summed. In the equations of Schacher 


and de Wette this introduces the following four terms into 


о увале 

giro eh 

[Ae ihe 

GR, Th?) 

CB, Th?) 
where h-h(.u, Az) j didla summing indices; and 
Г(а,х) «Ге ғар the incomplete gamma function. 

In order to develop a usable program based on Schacher 

and de Wette's equation it was necessary to obtain a simple 
subroutine for evaluating these functions repeatedly as h 
changed with A, апа 42. Unfortunately incomplete gamma 
functions are not widely used and the only subroutines 
available from other sources were those designed to tabulate 
Г (а, x) over a wide range of values of a and x. Since 


Г (а, х) is a transcedential whose "best! approximation is 
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an assymptotic series, tabulating programs are both 
cumbersome and slow. 

Three possibilities were considered for this program. 
First, adapt'enewot the existing? pregrams ferslimited use. 
This was tried but even after abstracting the desired 
parts from the best of these programs we were left with 
a double-precision program with an evaluation time, per h 
value, of about 0.5 sec. Since we contemplated hundreds 
of thousands of evaluations per lattice we ultimately 
discarded this approach. 

The second possibility considered was simply storing 
tabulated values over the range h-0-100 and interpolating 
as necessary. However to obtain program accuracies of 
i107? for S would have required storing 4000 or more 
values. ‘This, added to the large program storage already 
generated would have overflowed the internal storage of 
all but the largest computers available. Tape or disc 
storage would again increase run time too much so this 
approach was also discarded. 

The final possibility was adopted, that of formulating 
new approximations to the incomplete gamma function. 


Ultimately a set of power series approximations was developed, 
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using the Tau Method of emen to span the range 
fT h* -3-00. Another, long-known series, well convergent 
for small values of x, was used in the range Wh? =0-3 


The expansions are 


до (ПКМ 
M -R is KIM hoT О<хе3 


P (vex) 
,GY' E c, Q2" 34x 


where the C,;:, were generated using the Tau Method. The 
other functions, /(¥%,x) , etc were obtained from 7(%,x) 
by using'the recursion relation / (9+/,0) :0/7'(0,4) +092 . 
The Tau Method is well adopted to exploring several 
possibilities of approximating series on the computer. In 
this instance we examined a total of 240 possibilities, 
combining the four a values, -i 2, 8, i with 3 possible 
ranges of approximation, :<Х ,/<Х,.2<Ж and in each case 
examining series of 2, 3,  ..., 21 terms. Іп general, 
accuracy increased with increasing a value, increasing end 
point b, b« x, and increasing number of terms in the series. 


We used n vice Е ог $ because it matched the a value of 


SO C., Applied Analysis, 61-185, ‘Prentice 
Hall, Englewood Cliffs, N.J., 1961) 
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the expansion for O& x€3. A 17 term series was used 
since it exhausted machine accuracy of the computer used 


(+10711) 


ә апа any higher series would add no accuracy. 
Having settled on an approximation it was checked 

against a double precision routine known to be accurate to 

sup. The graph of the error of the approximation is 


shown in figure 13. 


7t 


2 А Г 
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Figure 13: 





Error curve for Incomplete Gamma F'unction 
approximation. sdbscissa is x, ordinate is 

* (%2,x)-1'(4a,x) where/’* is approximation and 

/" is actual value of Incomplete Gamma Function. 
Shaded area represents an oscillation of error 
of such a high frequency that it could not be 
faithfully represented on this scale. 
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APPENDED 


APPLICATION IROLBENZENEB 


During the preparation of this program R.M. Hexter 
of the Mellon Institute suggested an application to a long- 
standing problem in molecular spectroscopy. ‘This work has 
been completed and is detailed below as an example of applying 
the program to an actual problem. 

In calculating the electronic spectrum of benzene it 
is necessary to compute the transition moments of the free 
molecules. These moments depend onthe dipole-dipole sums 
of the crystal. Various workers computing these sums by 
various methods have obtained varying values for the sums. 
Craig and wale have set out the crystalline structure 
assumed and tabulated the values obtained by two sunmation 
procedures. ‘This program was used to check these values 
with results tabulated below. 

The benzene crystal model consists of an orthorhombic 
unit cell with four molecular orientations at the corners 
and faces. If we superimpose cartesian reference frames 


on each molecule we can identify our molecules as follows: 
9D. Craig and J. Walsh, J. Chem Soc 1958, 1613 


Е 


Molecular Unit Cell Crystal Direction Cosines 


Type Location Axis x y 2 

2 0,6,0 а .6365 -.2992 7109 
b .2000 „951 ‚2229 
© -.7450 0 . 6671 

II 5,0,7 а -.6365 .2992 -.7109 
b -.2000 -.9541 -.2229 
с =. 7450 0 16671 

III 0,4,4 a -. 6365 -29092 =.7109 
b .2000 .9541  =.2229 
C Т5 0; 0 Бб 

IV 2,1,0 а .6365 -.2992 27109 
b =. 2000  —J95h1 99229 


Where (a,b,c) forms the cartesian frame of the unit cell. 


The unit cell parameters are (at -3? C) 
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a4 7 7-460 A or =B=¥= 90° 
а; = 9.666 А 
аз = 7.035 А 


The problem then is to compute the interaction energy 
of a unit dipole in the i- direction in molecule I due to a 
lattice of unit dipoles in the j-direction in molecules of 
куре Ж, whére i |і = x, у, Zed — Ieee, III iv. 

In comparing figures obtained by this program with 
those tabulated by Craig and Walsh two difficulties arise. 
First dipole-dipole energies from this program work out in 


recipricol angstroms per coulomb squared while Craig and 


7% 


Walsh tabulate in recipricol centimeters (wave number) per 
electronic charge squared. This requires multiplying the 


results of this program by 1.161 X 10? 


to be comparable. 
Second Craig and Walsh assume an.infinite*sphéric x] 
crystal while this program assumes an infinite slab-shaped 


crystal. ‘Thus depolarization factors must be applied. The 


two systems are: 


Sperical Slab-shaped 
ATT 
£yx 7 3Vor £y, tO 
и 
Lyy ~ 3Vo Eyy tÒ 
E = 2.47: 
£j = Ezz Vol” 


Thus, to be comparable we must add = (477 /vol) = 


У 


"м ~. 


05017655170 ИЕ mu terms generated by this program and 
subtract - (477 /vol).= 0.0082585 to xx and xy terms 
generated by this program. (Where vol = volume of unit 
cell in A) 
GENERAL PRO&EDURE 

With these considerations in mind the computation 
becomes straight-forward if tedious. First we will compute 
electric field components at molecule I along the crystal 
axes a, b, c due to an assumed set of unit dipoles along 


the a, b, and c axes at molecular positions I, II, III, and 


T5 


% 
“. “р “= 


IV. This is done by calculating the field at the appropriate 
position of molecule I (given by XJI, XJ2, XJ3) within 
each of the I, II, III, and IV sublattice unit cells. Then 
we correct the depolarization factors and apply the conversion 
factor. Then, using the direction cosines of the type I 
molecule for each molecular axis x, y, and z (with respect 
to the a, b, c axes) we can get the x, y, and z components 
of field at I due to a, b, and c dipoles at all molecular 
positions. Then by multiplying by the direction cosines for 
each molecule respectively we get energies of the x, y, and 
z dipoles of molecule I due to x, y, and z dipoles at each 
molecular position. These we compare with similarly 
tabulated results of Craig and Walsh. 
DETAILS OF COMPUTATION 
To get the a, b, and c field components at molecule 


I we enter the program with: 


A1 = 7.460 XALPH = 90. 
A2 = 9.666 XBET = 90. 
A3 = 7.034 XGAM = 90. 

XDEL - 0. 


AEBX=REX=ABY=REY=AEZ=R ae 


N=20 


and for various molecules: 
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E IT III BA 
xai 0 п 0 4 
Gg 0 0 + 5 
X4J3 0 4 5 0 


Since the crystal is orthorhombic we get 570, ive 
The three nonzero outputs SUMXX, SUMYY, and SUM22 
here correspond to the a, b, and c:field components. ‘The 
computed S values, depolarization correction, conversion 
factor, and converted values are listed below: 


Fields at i-field 


I due to component 
molecules due to j- Computed Corrected 
of type dipoles Values Depolarization Values 
I aa 40102233 -.0082585 .0021648 
bb «OOM 720 -.0082585 .0070865 
ec .0115953 .0165170 .0049217 
TJ aa “02 LOZZO -.0082585 .0157635 
bb -.0187985 -.0082585 —.0270570 
се = „010052235 2096 1:950 40112935 
TIT aa -.0055470 -.0082585 “01238055 
bb .0231906 -.0082585 .0149321 
ес —-.0176436 40165170 =. 0011266 
puo aa 0080050 -.0082585 - 72072535 
bb 020918) -.0082585 .0126599 
ec — . 0289234 104651770 ~.012)06) 
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Molecule Axes 
I aa 
bb 
сс 
ІІ аа 
bb 
сс 
T IMP aa 
bb 
сс 
IV aa 
bb 
ce 


Corrected 
Value 


• обешена 
20070865 
0049217 


‚0157635 
ano Ob О 
.0112935 


‚0438055 
.01,9321 
. 0941966 


-00025935 
.0126599 
.0124064 


Now we need the x, y, 


Conversion 
Factor 


1. 161.302 


1.161х102 


Converted 
Values 


251.333 


=вег. lic 


571.409 


1830.142 
- 38 14817 
doped 75 


-1602.818 
17 ЗӘШЕШС 
-130.798 


1469.81} 
-1440.383 


and z components of these 


Meldesat тојевшент: І-І (хп), Шет (у, т) 2) (2,1), 


екс (т-а,Б,сф 


cesimes Гог weeenezewibthasespesibeto a, by and с} 


examples suffice: 


x - field components 


Fields at I 
due to dipoles Unit Call 
at Component 
" aa 
bb 
се 
ШІ ! 
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Here we need the molecule I direction 


Two 
Molecular 
Value Dir. Cos Component 
e DS .6365 т59. 97 
-822.'742 ‚2000 -164.5484 
571.409 =-.7450 -425.0997 


«y — field components 


Fields at I 
due to dipoles Unit Call Molecular 
at Component Value Dir Cos Component 
! 
! 
! 
Шев аа 1830.142 -.2992 =5 7.5976} 
bb ЕО ЗАР 1 -2997 . 1305 
ce 1311075 0 0 
Tq 


From these we can make up the energies by multiplying 
the molecular components by the direction cosines of the 
source dipole at the molecular position of interest. 

Since: we have tabulated the molecular components 
I-I (x, m) and I-II (y,m), if we want the energy of an 
x-dipole in molecule I due to a lattice of y-dipoles at 
position I, I-I (x,y), we need only multiply the molecular 
components I-I (x, m) by the y direction cosines of type 
I and sum. We show the computation of I-I (x,y) and 


I-II (y,z) below: 
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Field ja Direction Energy 


Component Axis Value Cosine Contribution 
Іу dipoles 
І-І (х, ) aa 159.973; -.2992 -47 . 8640 
bb -164.548 9541 -156.9956 
се =-425.6997 0 0 


I-I (x,y) energy -204.8596 


112 dipoles 


I-II (y, ) аа -547.57854 -.7109 289: 273 
bb-2997.1305 -.2229 668.0603 
сс 0 Сс”! 0 


I-II (y,z) energy 1057.3337 

These then are the numbers to be compared with 
those tabulated by Craig and Walsh. ‘Their table contains 
numbers obtained both from the Ewald-Kornfeld? method 
applied to an infinite sphere and from a direct molecule- 
molecule summation over a sphere of radius 20 A. Тһе 
comparison is shown below. Results (a) are from this 
program, (b) from the Ewald-Kornfeld method, and (c) 


from direct molecule-molecule summation. 


9 
Mitte) introduction to) Solid vacate ры есы, па кш 
347, (John Wiley and Sons, Inc., New York, 195) 
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Dipole 
Axes 


х, 


у, 


2; 


х, 


х, 


У, 


x 


Results 


(a) 
(b) 
(с) 


(а) 
(Ы) 
(с) 


(а) 
(Ы) 
(ке) 


(а) 
(5) 
(с) 


(а) 
(Ы) 
(с) 


(а) 
(Ы) 
(е) 


MET 


-386. 


-386 
-415 


(oon 


727 
790 


-240 
-341 
-375 


2a. 


205 
233 


ACC. 


207 
181. 


228 
225 
260 


06 


45 


e47 


86 


94 


43 


Molecule 


ІШІП 


-1110 


-112 
-92 


-2695. 


-2679 
-2629 


T 


TO 
20 


–9}7 
-947 
-934 


1339. 


1340 
1306 


-1057. 


-1058 
-1041 


n 


49 


‚96 


pe 


35 


Sites 


Т, СТО: 


-791. 


-791 
-769 


-1721 
-1723 
-1713 


-955. 


= 9 
=950 


-25. 


ES 
-8 


mio 


EU 
-818 


-2. 


s29 
-7 


29 


"OT 


10 


Б 


54 


76 


-728. 


-729 
-771 


1340 
ІРІ 
1405 


eue 


E 
-587 


Ou. 


275 
289 


79%. 


(99 
833 


306. 


S07 
323 


Ve 


рег 


РО 
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